!========================================================================================
! VERSION 1.00
!=========================================================================================
  module topology_m
!=========================================================================================
! module for topology variables

  implicit none
!-----------------------------------------------------------------------------------------
  include 'mpif.h'
!-----------------------------------------------------------------------------------------
! integers

  integer :: npes          !total number of processors
  integer :: myid          !rank of local processor
  integer :: xpes          !number of processors in x-direction
  integer :: xid           !id of local processor in x-direction
  integer :: ypes          !number of processors in y-direction
  integer :: yid           !id of local processor in y-direction
  integer :: zpes          !number of processors in z-direction
  integer :: zid           !id of local processor in z-direction

  integer :: term_status   !termination status (SDM)

  integer :: gcomm, xcomm, ycomm, zcomm

  integer :: neighbor(6)

  integer :: ierr, status(MPI_Status_size)
  integer :: hdf5_error


!----------------------------------------------------------------------------------------
  contains
!========================================================================================
  subroutine initialize_topology(io,npx,npy,npz)
!========================================================================================
! routine initializes some MPI stuff and the Cartesian MPI grid
!----------------------------------------------------------------------------------------
  implicit none

!----------------------------------------------------------------------------------------
! declarations passed in

  integer io
  integer nx, ny, nz, npx, npy, npz

! local declarations
  
  integer tmp
  integer mypx, mypy, mypz

  logical qflag

!----------------------------------------------------------------------------------------
! set term_status=0 for all processors

  term_status = 0
!----------------------------------------------------------------------------------------
!	initialize MPI stuff

  term_status=0

  call MPI_Initialized(qflag, ierr)

  if(.not.qflag) then
    call MPI_Init(ierr)
    call MPI_Comm_rank(MPI_COMM_WORLD, myid, ierr)
    call MPI_Comm_size(MPI_COMM_WORLD, npes, ierr)

!   create communicator duplicate for global calls

    call MPI_Comm_dup(MPI_COMM_WORLD, gcomm, ierr)

  endif

!----------------------------------------------------------------------------------------
! write header

  if(myid.eq.0) then
    write(io,*) 'Initializing topology module...'
    write(io,*)
  endif

!----------------------------------------------------------------------------------------
! check for npes compatibility

  if(npx*npy*npz.ne.npes) then
    if(myid.eq.0) then
      write(io,1) npx, npy, npz, npes
    endif
  endif

1 format(' npx*npy*npz is not equal to npes, npx = ',  &
            i5,' npy = ', i5, ' npz = ', i5, ' npes = ', i5)
!----------------------------------------------------------------------------------------
!	initialize Cartesian grid

  mypz = myid/(npx*npy)
  mypx = mod(myid-(mypz*npx*npy), npx)
  mypy = (myid-(mypz*npx*npy))/npx

  tmp = mypx-1 
  if(tmp.lt.0) then
    neighbor(1) = -1
  else 
    neighbor(1) = myid-1
  endif

  tmp = mypx+1 
  if (tmp.gt.npx-1) then
    neighbor(2) = -1
  else 
    neighbor(2) = myid+1
  endif

  tmp = mypy-1 
  if (tmp.lt.0) then
    neighbor(3) = -1
  else 
    neighbor(3) = myid-npx
  endif

  tmp = mypy+1 
  if (tmp.gt.npy-1) then
    neighbor(4) = -1
  else 
    neighbor(4) = myid+npx
  endif

  tmp = mypz-1 
  if (tmp.lt.0) then
    neighbor(5) = -1
  else 
    neighbor(5) = myid-(npx*npy)
  endif

  tmp = mypz+1 
  if (tmp.gt.npz-1) then
    neighbor(6) = -1
  else 
    neighbor(6) = myid+(npx*npy)
  endif

!----------------------------------------------------------------------------------------
!	create communicators for the x, y, and z directions

  call MPI_Comm_split(gcomm, mypy+1000*mypz, myid, xcomm,ierr)
  call MPI_Comm_split(gcomm, mypx+1000*mypz, myid, ycomm,ierr)
  call MPI_Comm_split(gcomm, mypx+1000*mypy, myid, zcomm,ierr)

  call MPI_Comm_rank(xcomm, xid, ierr)
  call MPI_Comm_size(xcomm, xpes,ierr)

  call MPI_Comm_rank(ycomm, yid, ierr)
  call MPI_Comm_size(ycomm, ypes,ierr)

  call MPI_Comm_rank(zcomm, zid, ierr)
  call MPI_Comm_size(zcomm, zpes,ierr)

  return
!----------------------------------------------------------------------------------------
  end subroutine initialize_topology

!-----------------------------------------------------------------------------------------
  end module topology_m